Analysis of Role of Dams on Recorded Dry Season Malaria Incidences in Kasungu
Clinton Nkolokosa
Last edited 11 June 2021
This geospatial statistical model uses routinely collected malaria case data, population data and remotely sensed data, such as open and vegetated water bodies, to estimate population living around open water bodies, expected malaria cases, and standardised morbidity ratio (SMR) of malaria. And ultimately, quantify the association between proximity to larval habitat and malaria risk in health facility catchment areas in Kasungu. The SMR compares the risk of morbidity in a population of interest with that of a standard population. In this case, our interest is to find out whether the number of dry season malaria cases in each catchment area are greater than we would expect given the malaria rate for the entire Kasungu district.
We do this by comparing what we observe (O) with what we would expect (E) if the risk of malaria was equal throughout Kasungu. The SMR statistical notation of catchment i can be written as follows: \[SMR_i = \frac{O_i}{E_i}\]
Buffers around waterbodies are created and then combined with population data in raster format to estimate the proprtion of catcment population living within 1km, 2km and 3km of water bodies. Subsequently, the observed malaria cases are modeled using Poisson regression to find out if living within various distances from water bodies is causing variability in malaria risk in Kasungu district. We hypothesize that the risk of being a case in a catchment is dependent on proximity to water bodies. The data used spans from 2017 to 2020 and was derived from digitized DHIS2 malaria records, accessibility mapping, aggregated population geospatial layer and TropWet tool in Google Earth Engine.
Load packages
Loading the R packages that will be used to read in, view, transform and model the malaria cases and spatial datasets.
library(spatialEco)
library(dplyr)
library(tidyverse)
library(ggplot2)
library(ggpubr)
library(plotly)
library(lubridate)
library(knitr)
library(raster)
library(rgdal)
library(rgeos)
library(sf)
library(sp)
library(tmap)
library(spdep)
library(maptools)
library(gridExtra)
library(grid)
library(exactextractr)
library(DataExplorer)
library(mapview)
`%>%` <- magrittr::`%>%`Tell R where the data is
here::here()## [1] "C:/Users/cnkolokosa/Documents/R/upscaled_2021_updated_May/upscaled_2021"Load datasets
The total dry season malaria cases recorded at health-care facilities in Kasungu from 2017 to 2019 are contained in the KasunguData.csv sourced from https://dhis2.health.gov.mw/. The kasungu_facility_catchments_2004.shp shapefile also contains the population and health information within each health-facility catchment area in Kasungu district.
The aggregated population raster layers for Malawi e.g.,ku_pop_2017_1km_aggregated.tif were downloaded from the Open Spatial and Demographic and Data Research website: https://www.worldpop.org/geodata/country?iso3=MWI. These layers estimate total number of people per grid-cell. The units are number of people per pixel with country totals adjusted to match the corresponding official United Nations population estimates. The datasets were downloaded in Geotiff at a resolution of 1km and are projected in Geographic Coordinate System, WGS84.
The kasungu_water.shpand water_bodies layers contain open and vegetated waterbodies polygons, detected using the Tropical Wetland Unmixing Tool (TropWet). TropWet is a Google Earth Engine hosted toolbox that uses the Landsat archive to map tropical wetlands and can be accessed through: https://www.aber.ac.uk/en/dges/research/earth-observation-laboratory/research/tropwet/
# 2017, 2018 and 2019 dry season malaria cases
dry_season_malaria_2015_2019 <- read.csv(here::here("data", "dry_season_malaria 2015-2019.csv"))
# Export 2017, 2018 and 2019 dry season malaria cases
write.csv(dry_season_malaria_2015_2019, file = "data/dry_season_malaria_2017_2019.csv")
write.table(dry_season_malaria_2015_2019,
file = "dry_season_malaria_2017_2019.txt",
sep = ",",
quote = FALSE,
row.names = FALSE)
# 2020 dry season malaria cases
ku_malaria_2020 <- read.csv(here::here("data", "dry_season_malaria_2020.csv"))
# Merge 2015 to 2019 dry season malaria case data with 2020 data
dry_season_malaria_2017_2020 <- cbind.data.frame(dry_season_malaria_2015_2019, ku_malaria_2020)
dry_season_malaria_2017_2020 <- dry_season_malaria_2017_2020[,c("FID", "Names", "dr_2017",
"dr_2018", "dr_2019", "dr_2020",
"LONGITU", "LATITUD")]
# Export 'dry season malaria 2017-2020' as csv
write.csv(dry_season_malaria_2017_2020, file = "data/dry_season_malaria_2017_2019.csv")
# Kasungu district boundary shapefile
kasungu_district <- sf::st_read(here::here("data", "kasungu_district.shp"))## Reading layer `kasungu_district' from data source `C:\Users\cnkolokosa\Documents\R\upscaled_2021_updated_May\upscaled_2021\data\kasungu_district.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 5 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 491272.7 ymin: 8494349 xmax: 609044.2 ymax: 8632164
## Projected CRS: WGS 84 / UTM zone 36S# Kasungu health facility catchments generated from accessibility mapping
malire_new <- sf::st_read(here::here("data", "new_catchments.shp")) %>%
sf::st_transform(32736) # reproject to WGS UTM Zone 36 South## Reading layer `new_catchments' from data source `C:\Users\cnkolokosa\Documents\R\upscaled_2021_updated_May\upscaled_2021\data\new_catchments.shp' using driver `ESRI Shapefile'
## Simple feature collection with 27 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 32.925 ymin: -13.61667 xmax: 34.00833 ymax: -12.375
## Geodetic CRS: WGS 84# Kasungu population raster layer
kasungu_population_2017 <- raster(here::here("data", "ku_pop_2017_1km_aggregated.tif"))
kasungu_population_2018 <- raster(here::here("data", "ku_pop_2018_1km_aggregated.tif"))
kasungu_population_2019 <- raster(here::here("data", "ku_pop_2019_1km_aggregated.tif"))
kasungu_population_2020 <- raster(here::here("data", "ku_pop_2020_1km_aggregated.tif"))
# Read in waterbodies polygons
dryseason_waterbodies_2017 <- sf::st_read(here::here("data", "water_bodies_2017.shp"))## Reading layer `water_bodies_2017' from data source `C:\Users\cnkolokosa\Documents\R\upscaled_2021_updated_May\upscaled_2021\data\water_bodies_2017.shp' using driver `ESRI Shapefile'
## Simple feature collection with 168 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 514497 ymin: 8495941 xmax: 603149.8 ymax: 8620169
## Projected CRS: WGS 84 / UTM zone 36Sdryseason_waterbodies_2018 <- sf::st_read(here::here("data", "kasungu_2018_water.shp"))## Reading layer `kasungu_2018_water' from data source `C:\Users\cnkolokosa\Documents\R\upscaled_2021_updated_May\upscaled_2021\data\kasungu_2018_water.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1105 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 496807.6 ymin: 8494693 xmax: 607913.8 ymax: 8607747
## Projected CRS: WGS 84 / UTM zone 36Sdryseason_waterbodies_2019 <- sf::st_read(here::here("data", "kasungu_2019_water.shp"))## Reading layer `kasungu_2019_water' from data source `C:\Users\cnkolokosa\Documents\R\upscaled_2021_updated_May\upscaled_2021\data\kasungu_2019_water.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1941 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 494197.2 ymin: 8494693 xmax: 607913.8 ymax: 8617573
## Projected CRS: WGS 84 / UTM zone 36Sdryseason_waterbodies_2020 <- sf::st_read(here::here("data", "water_bodies_2020.shp"))## Reading layer `water_bodies_2020' from data source `C:\Users\cnkolokosa\Documents\R\upscaled_2021_updated_May\upscaled_2021\data\water_bodies_2020.shp' using driver `ESRI Shapefile'
## Simple feature collection with 266 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 508985.6 ymin: 8495793 xmax: 585761.1 ymax: 8620169
## Projected CRS: WGS 84 / UTM zone 36S# Add a field ID to water bodies polygons
dryseason_waterbodies_2017$ID <- 1:nrow(dryseason_waterbodies_2017)
dryseason_waterbodies_2018$ID <- 1:nrow(dryseason_waterbodies_2018)
dryseason_waterbodies_2019$ID <- 1:nrow(dryseason_waterbodies_2019)
dryseason_waterbodies_2020$ID <- 1:nrow(dryseason_waterbodies_2020)View the dry season malaria case data
We observe that Kasungu district has 30 health facilities classified as dispensary, health centre, district hospital and rural hospital, and the highest malaria cases were recorded at Kasungu District Hospital.
dry_season_malaria_2017_2020 %>%
summary()## FID Names dr_2017 dr_2018
## Min. : 1.00 Length:36 Min. : 427 Min. : 749
## 1st Qu.: 9.75 Class :character 1st Qu.: 2196 1st Qu.: 2424
## Median :18.50 Mode :character Median : 3132 Median : 3357
## Mean :18.50 Mean : 3586 Mean : 3909
## 3rd Qu.:27.25 3rd Qu.: 3993 3rd Qu.: 4684
## Max. :36.00 Max. :16289 Max. :15821
## NA's :6 NA's :6
## dr_2019 dr_2020 LONGITU LATITUD
## Min. : 533 Min. : 0 Min. :33.18 Min. :-13.57
## 1st Qu.: 1748 1st Qu.: 1650 1st Qu.:33.38 1st Qu.:-13.25
## Median : 2556 Median : 2698 Median :33.50 Median :-12.98
## Mean : 2880 Mean : 3404 Mean :33.52 Mean :-12.99
## 3rd Qu.: 3284 3rd Qu.: 4722 3rd Qu.:33.68 3rd Qu.:-12.79
## Max. :10721 Max. :16435 Max. :33.87 Max. :-12.42
## NA's :6 NA's :6 NA's :6dry_season_malaria_2017_2020 %>%
plotly::plot_ly(y = ~Names,
x = ~dr_2017,
type = "bar",
orientation = 'h',
name = "2017") %>%
plotly::add_trace(x = ~ dr_2018,
name = "2018") %>%
plotly::add_trace(x = ~ dr_2019,
name = "2019") %>%
plotly::add_trace(x = ~ dr_2020,
name = "2020") %>%
plotly::layout(xaxis = list(title = "Total malaria cases"),
yaxis = list(title = " "),
hovermode = "compare",
margin = list(b = 10,
t = 10,
pad = 2))Fig.1 The total malaria cases recorded at each health-care facility in Kasungu district
Kasungu health-care facilities and their catchment areas
Heath facility catchment area is the area from which a health facility attracts patients. The new health facility catchments polygon was generated from generic accessibility mapping script adapted from https://malariaatlas.org/wp-content/uploads/accessibility/R_generic_accessibilty_mapping_script.r The script requires two user supplied datasets: the 2015 friction surface, which is available here: http://www.map.ox.ac.uk/accessibility_to_cities/, and a user-supplied .csv of points dry_season_malaria_2017_2020. The accumulated cost algorithm accCost and r.Cost algorithm in GRASS GIS were run to make the final output map of new health facility catchment boundaries.
# Using the complete.cases() function to select health centres with complete
# longitude and latitude coordinates.
zipatala_aggregated <- dry_season_malaria_2017_2020[complete.cases(dry_season_malaria_2017_2020),]
# Aggregate health facilities close to each other:
# a) Kasalika Health Centre and Kasungu District Hospital,
# b) Bua and Mziza Health Centres, and
# c) Kaluluma and Nkhamenya Rural Hospitals in order to
# generate catchment areas that are geographically correct
zipatala_aggregated$dr_2017[which(
zipatala_aggregated$Names == "Kasungu District Hospital")] <- zipatala_aggregated$dr_2017[which(
zipatala_aggregated$Names == "Kasungu District Hospital")] + zipatala_aggregated$dr_2017[which(
zipatala_aggregated$Names == "Kasalika Health Centre")]
zipatala_aggregated$dr_2018[which(
zipatala_aggregated$Names == "Kasungu District Hospital")] <- zipatala_aggregated$dr_2018[which(
zipatala_aggregated$Names == "Kasungu District Hospital")] + zipatala_aggregated$dr_2018[which(
zipatala_aggregated$Names == "Kasalika Health Centre")]
zipatala_aggregated$dr_2019[which(
zipatala_aggregated$Names == "Kasungu District Hospital")] <- zipatala_aggregated$dr_2019[which(
zipatala_aggregated$Names == "Kasungu District Hospital")] + zipatala_aggregated$dr_2019[which(
zipatala_aggregated$Names == "Kasalika Health Centre")]
zipatala_aggregated$dr_2020[which(
zipatala_aggregated$Names == "Kasungu District Hospital")] <- zipatala_aggregated$dr_2020[which(
zipatala_aggregated$Names == "Kasungu District Hospital")] + zipatala_aggregated$dr_2020[which(
zipatala_aggregated$Names == "Kasalika Health Centre")]
zipatala_aggregated$dr_2017[which(
zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] <- zipatala_aggregated$dr_2017[which(
zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] + zipatala_aggregated$dr_2017[which(
zipatala_aggregated$Names == "Kaluluma Rural Hospital")]
zipatala_aggregated$dr_2018[which(
zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] <- zipatala_aggregated$dr_2018[which(
zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] +zipatala_aggregated$dr_2018[which(
zipatala_aggregated$Names == "Kaluluma Rural Hospital")]
zipatala_aggregated$dr_2019[which(
zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] <- zipatala_aggregated$dr_2019[which(
zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] + zipatala_aggregated$dr_2019[which(
zipatala_aggregated$Names == "Kaluluma Rural Hospital")]
zipatala_aggregated$dr_2020[which(
zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] <- zipatala_aggregated$dr_2020[which(
zipatala_aggregated$Names == "Nkhamenya Rural Hospital")] + zipatala_aggregated$dr_2020[which(
zipatala_aggregated$Names == "Kaluluma Rural Hospital")]
zipatala_aggregated$dr_2017[which(
zipatala_aggregated$Names == "Mziza Health Centre")] <- zipatala_aggregated$dr_2017[which(
zipatala_aggregated$Names == "Mziza Health Centre")] + zipatala_aggregated$dr_2017[which(
zipatala_aggregated$Names == "Bua Health Centre")]
zipatala_aggregated$dr_2018[which(
zipatala_aggregated$Names == "Mziza Health Centre")] <- zipatala_aggregated$dr_2017[which(
zipatala_aggregated$Names == "Mziza Health Centre")] + zipatala_aggregated$dr_2017[which(
zipatala_aggregated$Names == "Bua Health Centre")]
zipatala_aggregated$dr_2019[which(
zipatala_aggregated$Names == "Mziza Health Centre")] <- zipatala_aggregated$dr_2017[which(
zipatala_aggregated$Names == "Mziza Health Centre")] + zipatala_aggregated$dr_2017[which(
zipatala_aggregated$Names == "Bua Health Centre")]
zipatala_aggregated$dr_2020[which(
zipatala_aggregated$Names == "Mziza Health Centre")] <- zipatala_aggregated$dr_2017[which(
zipatala_aggregated$Names == "Mziza Health Centre")] + zipatala_aggregated$dr_2017[which(
zipatala_aggregated$Names == "Bua Health Centre")]
# Drop out the other health facilities
zipatala_aggregated <- zipatala_aggregated %>%
dplyr::filter(Names != "Kasalika Health Centre",
Names != "Bua Health Centre",
Names != "Kaluluma Rural Hospital")
# write.csv(zipatala_aggregated, "data/health_facilities_aggregated.csv")
# Convert to csv spatial to spatial object
health_facility_aggr_sf <- sf::st_as_sf(zipatala_aggregated,
coords = c("LONGITU", "LATITUD"),
crs = 4326, agr = "constant")
# st_write(health_facility_aggr_sf, "data/health_facilities_aggregated.shp")View location of the health facilities in the new catchment areas
# Plot map
tm_shape(malire_new)+
tm_polygons()+
tm_shape(health_facility_aggr_sf)+
tm_dots(size = .3,
col = "blue",
alpha = 0.5)+
tm_text("Names",
size = .3,
just = "top",
col = "black",
remove.overlap = TRUE)+
tm_layout(frame = FALSE,
title = "New Kasungu health facility \n catchment boundaries",
title.size = .8,
title.position = c("left", "top"))+
tm_compass(position=c("right", "top"))+
tm_scale_bar(breaks = c(0, 10, 20),
text.size = .5)
Fig 2. Kasungu health-care facilities and catchment areas
Kasungu district estimated population per grid-cell
# Take a glimpse at the WorldPop raster layers
kasungu_population_2017## class : RasterLayer
## dimensions : 150, 128, 19200 (nrow, ncol, ncell)
## resolution : 920.0898, 918.7667 (x, y)
## extent : 491272.7, 609044.2, 8494349, 8632164 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs
## source : C:/Users/cnkolokosa/Documents/R/upscaled_2021_updated_May/upscaled_2021/data/ku_pop_2017_1km_aggregated.tif
## names : ku_pop_2017_1km_aggregatedkasungu_population_2018## class : RasterLayer
## dimensions : 150, 128, 19200 (nrow, ncol, ncell)
## resolution : 920.0898, 918.7667 (x, y)
## extent : 491272.7, 609044.2, 8494349, 8632164 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs
## source : C:/Users/cnkolokosa/Documents/R/upscaled_2021_updated_May/upscaled_2021/data/ku_pop_2018_1km_aggregated.tif
## names : ku_pop_2018_1km_aggregated
## values : 0, 6253.557 (min, max)kasungu_population_2019## class : RasterLayer
## dimensions : 150, 128, 19200 (nrow, ncol, ncell)
## resolution : 920.0898, 918.7667 (x, y)
## extent : 491272.7, 609044.2, 8494349, 8632164 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs
## source : C:/Users/cnkolokosa/Documents/R/upscaled_2021_updated_May/upscaled_2021/data/ku_pop_2019_1km_aggregated.tif
## names : ku_pop_2019_1km_aggregated
## values : 0, 6483.727 (min, max)kasungu_population_2020## class : RasterLayer
## dimensions : 150, 128, 19200 (nrow, ncol, ncell)
## resolution : 920.0898, 918.7667 (x, y)
## extent : 491272.7, 609044.2, 8494349, 8632164 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs
## source : C:/Users/cnkolokosa/Documents/R/upscaled_2021_updated_May/upscaled_2021/data/ku_pop_2020_1km_aggregated.tif
## names : ku_pop_2020_1km_aggregated
## values : 0, 7949.033 (min, max)# Helper function to create a raster population map
create.population.map <- function(population.raster, title){
# raster population map
# arguments:
# population.raster: aggregated population raster layer from WorldPop
# legend.title: legend title
# returns:
# a tmap-element (plots a map)
tm_shape(population.raster)+
tm_raster(palette = "YlOrBr",
title = title,
breaks = c(0,100,200,400,600,800,1000,2000,4000,6000,8000))+
tm_layout(legend.position = c("right", "bottom"),
frame = FALSE)+
tm_scale_bar(position = c("left", "bottom"))
}
# Set to static map
tmap_mode("plot")
estimated_pop_2017 <- create.population.map(kasungu_population_2017, title = "2017 Population")
estimated_pop_2018 <- create.population.map(kasungu_population_2018, title = "2018 Population")
estimated_pop_2019 <- create.population.map(kasungu_population_2019, title = "2019 Population")
estimated_pop_2020 <- create.population.map(kasungu_population_2020, title = "2020 Population")
# Layout the maps
tmap_arrange(estimated_pop_2017, estimated_pop_2018, estimated_pop_2019, estimated_pop_2020, nrow = 2) 
Fig.3 Estimated total number of people per 1km grid-cell
Assign dry season malaria cases and population density to new health facility catchments
The WorldPop aggregated population e.g. kasungu_population_2017.tif, and DHIS2 malaria dry_season_malaria_2017_2020 datasets are assigned to the new health facility catchments.
# Helper function that assigns malaria data from health facilities to their catchments areas ----------------
assign.malaria.data <- function(catchment_boundary, malaria_data){
# arguments:
# catchment_boundary: sf polygon object of new catchment boundaries
# malaria_data: sf point object with a data frame containing the dry season malaria cases
# returns:
# catchments_malaria_sf: sf polygon object with a data frame containing dry season malaria cases
# Convert sf objects to spatial
catchment_shp <- as(catchment_boundary, "Spatial")
malaria_shp <- as(malaria_data, "Spatial")
# Match CRS
malaria_shp <- spTransform(malaria_shp, crs(catchment_shp))
# Overlay aggregated health facility points and extract 2017 - 2020 malaria cases
# Using 'point.in.poly' to return a point spatial object, in this case location of health facilities
# and estimated population instead of sp::over function, which simply returns
# a data frame, with the same no. rows.
# Argument 'sp = TRUE' returns an sp class object, else returns sf class object
# Joining the malaria and population dataset using only 'merge' function can't work due to
# non-unique columns and differences in row numbers
hospitals_in_catchment <- spatialEco::point.in.poly(malaria_shp, catchment_shp, sp = TRUE)
# Add the extracted ID, health facility names and dry season malaria cases to
# the health facility catchments (hfc)
hfc_malaria_shp <- merge(catchment_shp, hospitals_in_catchment, by.x = "DN", by.y = "FID")
# Convert the shapefile containing malaria data to sf-object
hfc_malaria_sf <- sf::st_as_sf(hfc_malaria_shp)
# Tidy the data by dropping columns not needed
catchment_malaria <- hfc_malaria_sf %>%
dplyr::select(-c(coords.x1, coords.x2))
return(out = catchment_malaria)
}
# Invoking the function ----------------------------------------------------------------------------------
malaria_by_catchment <- assign.malaria.data(malire_new, health_facility_aggr_sf)Assign population data to the health catchment areas
# Helper unction to extract population from WorldPop raster file and assign ---------------------------
# the values to the new catchments.
extract.pop.values <- function(kasungu_pop_raster, catchments){
# function to extract population from raster file and assign the population to catchments
# arguments:
# kasungu_pop_raster: population raster file clipped to Kasungu district
# catchments: shapefile containing the polygons that we wish to use as boundaries
# returns:
# catchments_malaria_pop_sf: sf polygon object containing malaria and population data
# convert from sf to sp
catchments_sp <- as(catchments, "Spatial")
# Match extent i.e projection
catchments_sp <- spTransform(catchments_sp, proj4string(kasungu_pop_raster))
# Crop and mask the population raster to exclude Kasungu National Park
pop_raster_clip <- raster::mask(raster::crop(kasungu_pop_raster, extent(catchments_sp)), catchments_sp)
# Extracting zonal statistics from a population raster layer.
# The population raster is a continuous gridded surface layer that has an
# estimated population density value to every square in their grid.
# The population values are then summed and apportioned to the catchment polygons
# catchments_malaria_pop <- catchments %>%
# dplyr::mutate(pop = round(raster::extract(pop_raster_clip, catchments, fun = sum, na.rm = TRUE)))
pop_by_catchment <- round(raster::extract(pop_raster_clip, catchments, fun = sum, na.rm = TRUE))
pop_by_catchment_df <- pop_by_catchment %>%
# apply unlist to the lists to have vectors as the list elements
lapply(unlist) %>%
# convert vectors to data.frames
lapply(as_tibble) %>%
# combine the list of data.frames
bind_rows(., .id = "rowID") %>%
# rename the value variable
dplyr::rename(pop = value)
# Add row ID to column to catchment layer
catchments$rowID <- 1:nrow(catchments)
# Merge catchment areas and population data
pop_by_catchments <- merge(catchments, pop_by_catchment_df, by.x = "rowID", by.y = "rowID")
# Cleaning 'Inf' values
pop_by_catchments %>%
dplyr::mutate_if(is.numeric, list(~na_if(., Inf))) %>%
dplyr::mutate_if(is.numeric, list(~na_if(., -Inf)))
return(out = pop_by_catchments)
}
# Invoking the function ---------------------------------------------------------------------------------------
malaria_pop_by_catchment_2017 <- extract.pop.values(kasungu_population_2017, malaria_by_catchment)
malaria_pop_by_catchment_2018 <- extract.pop.values(kasungu_population_2018, malaria_by_catchment)
malaria_pop_by_catchment_2019 <- extract.pop.values(kasungu_population_2019, malaria_by_catchment)
malaria_pop_by_catchment_2020 <- extract.pop.values(kasungu_population_2020, malaria_by_catchment)View Kasungu population by catchment maps
Estimated total number of people within health facility catchment areas.
# Helper function to create maps of estimated population by catchment areas --------------------------------
create.population.map <- function(catchment.area,
variable = "pop",
title,
legend.title = "Estimated \n population"){
# estimated population map
# catchment.area: estimated population layer from nachulu function
# variable: variable name (as character, in qoutes)
# title: map title in quotes
# legend.title: legend title in qoutes
# returns:
# a tmap-element (plots a map)
tm_shape(catchment.area)+
tm_fill(col = variable,
breaks = c(0, 13000, 19000, 27000, 35000, 70000, 140000, 200000),
palette = "YlOrBr",
title = legend.title)+
tm_borders(col = "grey",
lwd = 0.4)+
tm_layout(legend.position = c(0.75, "bottom"),
legend.text.size = 0.6,
legend.title.size = 0.8,
frame = FALSE)+
tm_credits(title,
position = c(0.3, 0.8),
size = 1)
}
# Invoking the function --------------------------------------------------------------------------------
pop_by_catchment_2017 <- create.population.map(malaria_pop_by_catchment_2017, title = "2017")
pop_by_catchment_2018 <- create.population.map(malaria_pop_by_catchment_2018, title = "2018")
pop_by_catchment_2019 <- create.population.map(malaria_pop_by_catchment_2019, title = "2019")
pop_by_catchment_2020 <- create.population.map(malaria_pop_by_catchment_2020, title = "2020")
tmap::tmap_arrange(pop_by_catchment_2017, pop_by_catchment_2018,
pop_by_catchment_2019, pop_by_catchment_2020, ncol = 2)
Fig. 4: Estimated population by health facility catchment areas
Population density by catchment
# Helper function to calculate population density by catchment -----------------
calculate.population.density <- function(pop.data){
# Convert to spatial object
pop.sp <- as(pop.data, "Spatial")
# Calculate area of catchment polygon in square kilometres
pop.sp$area_sqkm <- round(rgeos::gArea(pop.sp, byid = TRUE) / (1000 * 1000))
# Calculate population density
pop.sp$pop_density <- round(pop.sp$pop / pop.sp$area_sqkm)
# Convert back to sf object
pop.sf <- sf::st_as_sf(pop.sp)
return(pop.sf)
}
# Invoking function ------------------------------------------------------------
pop_density_2017 <- calculate.population.density(malaria_pop_by_catchment_2017)
pop_density_2018 <- calculate.population.density(malaria_pop_by_catchment_2018)
pop_density_2019 <- calculate.population.density(malaria_pop_by_catchment_2019)
pop_density_2020 <- calculate.population.density(malaria_pop_by_catchment_2020)
# Helper function to create population density maps ----------------------------
create.pop.density.map <- function(pop.density.data,
variable = "pop_density",
title = NA,
legend.title = "Population \ndensity/sqkm"){
tm_shape(pop.density.data)+
tm_fill(col = variable,
breaks = c(0, 50, 100, 150, 200, 250, 300, 350),
palette = "-magma",
title = legend.title)+
tm_borders()+
tm_layout(legend.position = c(0.75, "bottom"),
legend.text.size = 0.6,
legend.title.size = 0.8,
frame = FALSE)+
tm_credits(title,
position = c(0.3, 0.8),
size = 1)
}
# Invoking function ------------------------------------------------------------
pop_density_2017_map <- create.pop.density.map(pop_density_2017, title = "2017")
pop_density_2018_map <- create.pop.density.map(pop_density_2018, title = "2018")
pop_density_2019_map <- create.pop.density.map(pop_density_2019, title = "2019")
pop_density_2020_map <- create.pop.density.map(pop_density_2020, title = "2020")
# Layout maps
tmap::tmap_arrange(pop_density_2017_map, pop_density_2018_map,
pop_density_2019_map, pop_density_2020_map, ncol = 2)
Fig. 5: Estimated population density by health facility catchment areas
Calculate the expected number of cases for each catchment area
The expected number of dry season malaria cases in catchment i are calculated as the observed risk (r) of malaria i.e. the total number of malaria cases in Kasungu district divided by the total population of the district, multiplied by the number of people in the catchment area: \[E_i = \frac{\sum_i O_i}{\sum_i N_i}\times N_i\]
The expected number of dry season malaria cases are calculated under the assumption that there is no spatial variation in risk, i.e., no difference in infection rates between the catchment areas.
# Calculate expected malaria cases --------------------------------------------------------------
expected_malaria_2017 <- malaria_pop_by_catchment_2017 %>%
dplyr::rename(
observed_2017 = dr_2017,
pop_2017 = pop) %>%
dplyr::mutate(
expected_2017 = round(sum(observed_2017)/sum(pop_2017, na.rm = TRUE)*pop_2017))
expected_malaria_2018 <- malaria_pop_by_catchment_2018 %>%
dplyr::rename(
observed_2018 = dr_2018,
pop_2018 = pop) %>%
dplyr::mutate(
expected_2018 = round(sum(observed_2018)/sum(pop_2018, na.rm = TRUE)*pop_2018))
expected_malaria_2019 <- malaria_pop_by_catchment_2019 %>%
dplyr::rename(
observed_2019 = dr_2019,
pop_2019 = pop) %>%
dplyr::mutate(
expected_2019 = round(sum(observed_2019)/sum(pop_2019, na.rm = TRUE)*pop_2019))
expected_malaria_2020 <- malaria_pop_by_catchment_2020 %>%
dplyr::rename(
observed_2020 = dr_2020,
pop_2020 = pop) %>%
dplyr::mutate(
expected_2020 = round(sum(observed_2020)/sum(pop_2020, na.rm = TRUE)*pop_2020))Calculate the Standardised Morbidity Ratio of malaria incidences for each catchment area
The SMR compares the risk of morbidity in a population of interest with that of a standard population. In this case, our interest is to find out whether the number of dry season malaria cases in each catchment area are greater than we would expect given the malaria rate for the entire Kasungu district.
We do this by comparing what we observe (O) with what we would expect (E) if the risk of malaria was equal throughout Kasungu. The SMR statistical notation of catchment i can be written as follows: \[SMR_i = \frac{O_i}{E_i}\]
# Calculate Standardised Morbidity Ratio (SMR) -------------------------------------
SMR_2017 <- expected_malaria_2017 %>%
dplyr::mutate(SMR = round(observed_2017/expected_2017, 1)) %>%
dplyr::select(rowID,Names, pop_2017, observed_2017, expected_2017, SMR)
SMR_2018 <- expected_malaria_2018 %>%
dplyr::mutate(SMR = round(observed_2018/expected_2018, 1)) %>%
dplyr::select(rowID, Names, pop_2018, observed_2018, expected_2018, SMR)
SMR_2019 <- expected_malaria_2019 %>%
dplyr::mutate(SMR = round(observed_2019/expected_2019, 1)) %>%
dplyr::select(rowID, Names, pop_2019, observed_2019, expected_2019, SMR)
SMR_2020 <- expected_malaria_2020 %>%
dplyr::mutate(SMR = round(observed_2020/expected_2020, 1)) %>%
dplyr::select(rowID, Names, pop_2020, observed_2020, expected_2020, SMR)
# Create SMR tables ------------------------------------------------------------
SMR_table_2017 <- SMR_2017 %>%
dplyr::as_tibble() %>%
dplyr::select(-rowID, -geometry) %>%
kable %>%
kableExtra::kable_styling(full_width = FALSE)
SMR_table_2018 <- SMR_2018 %>%
dplyr::as_tibble() %>%
dplyr::select(-rowID, -geometry) %>%
kable %>%
kableExtra::kable_styling(full_width = FALSE)
SMR_table_2019 <- SMR_2019 %>%
dplyr::as_tibble() %>%
dplyr::select(-rowID, -geometry) %>%
kable %>%
kableExtra::kable_styling(full_width = FALSE)
SMR_table_2020 <- SMR_2020 %>%
dplyr::as_tibble() %>%
dplyr::select(-rowID, -geometry) %>%
kable %>%
kableExtra::kable_styling(full_width = FALSE)
SMR_table_2017| Names | pop_2017 | observed_2017 | expected_2017 | SMR |
|---|---|---|---|---|
| Lodjwa Health Centre | 9923 | 1161 | 1372 | 0.8 |
| Nkhamenya Rural Hospital | 40154 | 3639 | 5553 | 0.7 |
| Newa Mpasazi Health Centre | 13879 | 427 | 1919 | 0.2 |
| Mpepa /Chisinga Health Centre | 27459 | 2784 | 3797 | 0.7 |
| Mnyanja Health Centre | 39950 | 3298 | 5524 | 0.6 |
| Simlemba Health Centre | 26999 | 2197 | 3734 | 0.6 |
| Ofesi Health Centre | 28098 | 4036 | 3886 | 1.0 |
| Chulu Health Centre | 27906 | 5801 | 3859 | 1.5 |
| Kapelula Health Centre | 35727 | 4138 | 4940 | 0.8 |
| Livwezi Health Centre | 22009 | 1307 | 3043 | 0.4 |
| Gogode Dispensary | 13061 | 2745 | 1806 | 1.5 |
| Dwangwa Dispensary | 32704 | 2192 | 4522 | 0.5 |
| Chamama Health Facility | 20026 | 1878 | 2769 | 0.7 |
| Wimbe Health Centre | 11864 | 5660 | 1641 | 3.4 |
| Chinyama | 12768 | 2292 | 1766 | 1.3 |
| Mdunga Health Centre | 18177 | 2045 | 2514 | 0.8 |
| Mtunthama Health Centre | 18744 | 3622 | 2592 | 1.4 |
| Kasungu District Hospital | 143490 | 20817 | 19842 | 1.0 |
| Chamwabvi Health Centre | 35353 | 3601 | 4889 | 0.7 |
| Linyangwa Health Centre | 17772 | 3359 | 2458 | 1.4 |
| Kawamba Health Centre | 22865 | 5808 | 3162 | 1.8 |
| Mziza Health Centre | 44189 | 7352 | 6111 | 1.2 |
| Kamboni Health Centre | 21226 | 3824 | 2935 | 1.3 |
| Khola Health Centre | 16956 | 2195 | 2345 | 0.9 |
| Santhe Health Centre | 6096 | 5838 | 843 | 6.9 |
| Anchor Farm | 48861 | 2966 | 6757 | 0.4 |
| Mkhota Health Centre | 21621 | 2586 | 2990 | 0.9 |
SMR_table_2018| Names | pop_2018 | observed_2018 | expected_2018 | SMR |
|---|---|---|---|---|
| Lodjwa Health Centre | 10281 | 1151 | 1588 | 0.7 |
| Nkhamenya Rural Hospital | 41642 | 4540 | 6430 | 0.7 |
| Newa Mpasazi Health Centre | 14248 | 749 | 2200 | 0.3 |
| Mpepa /Chisinga Health Centre | 28488 | 3602 | 4399 | 0.8 |
| Mnyanja Health Centre | 41856 | 2864 | 6464 | 0.4 |
| Simlemba Health Centre | 27455 | 2515 | 4240 | 0.6 |
| Ofesi Health Centre | 29002 | 3446 | 4479 | 0.8 |
| Chulu Health Centre | 28832 | 4502 | 4452 | 1.0 |
| Kapelula Health Centre | 37630 | 5338 | 5811 | 0.9 |
| Livwezi Health Centre | 22544 | 2361 | 3481 | 0.7 |
| Gogode Dispensary | 13368 | 4138 | 2064 | 2.0 |
| Dwangwa Dispensary | 33534 | 2394 | 5178 | 0.5 |
| Chamama Health Facility | 20372 | 1750 | 3146 | 0.6 |
| Wimbe Health Centre | 11814 | 5010 | 1824 | 2.7 |
| Chinyama | 13138 | 2116 | 2029 | 1.0 |
| Mdunga Health Centre | 18928 | 2923 | 2923 | 1.0 |
| Mtunthama Health Centre | 19074 | 5308 | 2945 | 1.8 |
| Kasungu District Hospital | 147175 | 20314 | 22727 | 0.9 |
| Chamwabvi Health Centre | 36167 | 4027 | 5585 | 0.7 |
| Linyangwa Health Centre | 18032 | 3268 | 2785 | 1.2 |
| Kawamba Health Centre | 22902 | 6462 | 3537 | 1.8 |
| Mziza Health Centre | 46208 | 10841 | 7136 | 1.5 |
| Kamboni Health Centre | 21430 | 4745 | 3309 | 1.4 |
| Khola Health Centre | 17315 | 2802 | 2674 | 1.0 |
| Santhe Health Centre | 6244 | 7267 | 964 | 7.5 |
| Anchor Farm | 49871 | 3142 | 7701 | 0.4 |
| Mkhota Health Centre | 22167 | 5920 | 3423 | 1.7 |
SMR_table_2019| Names | pop_2019 | observed_2019 | expected_2019 | SMR |
|---|---|---|---|---|
| Lodjwa Health Centre | 10608 | 1713 | 1200 | 1.4 |
| Nkhamenya Rural Hospital | 43293 | 4537 | 4896 | 0.9 |
| Newa Mpasazi Health Centre | 14780 | 809 | 1672 | 0.5 |
| Mpepa /Chisinga Health Centre | 29456 | 4248 | 3331 | 1.3 |
| Mnyanja Health Centre | 43783 | 3148 | 4952 | 0.6 |
| Simlemba Health Centre | 28076 | 2339 | 3175 | 0.7 |
| Ofesi Health Centre | 30065 | 4079 | 3400 | 1.2 |
| Chulu Health Centre | 29731 | 5470 | 3362 | 1.6 |
| Kapelula Health Centre | 39747 | 2558 | 4495 | 0.6 |
| Livwezi Health Centre | 22945 | 787 | 2595 | 0.3 |
| Gogode Dispensary | 13641 | 2200 | 1543 | 1.4 |
| Dwangwa Dispensary | 34415 | 2553 | 3892 | 0.7 |
| Chamama Health Facility | 20701 | 1508 | 2341 | 0.6 |
| Wimbe Health Centre | 11855 | 3041 | 1341 | 2.3 |
| Chinyama | 13475 | 1770 | 1524 | 1.2 |
| Mdunga Health Centre | 19960 | 1008 | 2257 | 0.4 |
| Mtunthama Health Centre | 19385 | 2763 | 2192 | 1.3 |
| Kasungu District Hospital | 151079 | 13450 | 17086 | 0.8 |
| Chamwabvi Health Centre | 36899 | 1686 | 4173 | 0.4 |
| Linyangwa Health Centre | 18279 | 3330 | 2067 | 1.6 |
| Kawamba Health Centre | 23041 | 4402 | 2606 | 1.7 |
| Mziza Health Centre | 48340 | 10841 | 5467 | 2.0 |
| Kamboni Health Centre | 21509 | 2770 | 2433 | 1.1 |
| Khola Health Centre | 17761 | 2708 | 2009 | 1.3 |
| Santhe Health Centre | 6435 | 5210 | 728 | 7.2 |
| Anchor Farm | 50995 | 1978 | 5767 | 0.3 |
| Mkhota Health Centre | 22677 | 2162 | 2565 | 0.8 |
SMR_table_2020| Names | pop_2020 | observed_2020 | expected_2020 | SMR |
|---|---|---|---|---|
| Lodjwa Health Centre | 13081 | 1788 | 1603 | 1.1 |
| Nkhamenya Rural Hospital | 53692 | 8539 | 6581 | 1.3 |
| Newa Mpasazi Health Centre | 18311 | 2182 | 2244 | 1.0 |
| Mpepa /Chisinga Health Centre | 36317 | 5186 | 4451 | 1.2 |
| Mnyanja Health Centre | 54649 | 6117 | 6698 | 0.9 |
| Simlemba Health Centre | 34240 | 5310 | 4197 | 1.3 |
| Ofesi Health Centre | 37240 | 2323 | 4564 | 0.5 |
| Chulu Health Centre | 36638 | 7160 | 4491 | 1.6 |
| Kapelula Health Centre | 50214 | 7297 | 6154 | 1.2 |
| Livwezi Health Centre | 27786 | 1028 | 3406 | 0.3 |
| Gogode Dispensary | 16681 | 2767 | 2045 | 1.4 |
| Dwangwa Dispensary | 42282 | 2869 | 5182 | 0.6 |
| Chamama Health Facility | 25248 | 635 | 3095 | 0.2 |
| Wimbe Health Centre | 14367 | 2233 | 1761 | 1.3 |
| Chinyama | 16463 | 1605 | 2018 | 0.8 |
| Mdunga Health Centre | 25108 | 3169 | 3077 | 1.0 |
| Mtunthama Health Centre | 23501 | 1882 | 2880 | 0.7 |
| Kasungu District Hospital | 185282 | 19393 | 22709 | 0.9 |
| Chamwabvi Health Centre | 45106 | 1128 | 5528 | 0.2 |
| Linyangwa Health Centre | 22144 | 4380 | 2714 | 1.6 |
| Kawamba Health Centre | 27961 | 7073 | 3427 | 2.1 |
| Mziza Health Centre | 60510 | 10841 | 7416 | 1.5 |
| Kamboni Health Centre | 25750 | 4665 | 3156 | 1.5 |
| Khola Health Centre | 21929 | 3426 | 2688 | 1.3 |
| Santhe Health Centre | 7917 | 4891 | 970 | 5.0 |
| Anchor Farm | 62633 | 1665 | 7677 | 0.2 |
| Mkhota Health Centre | 27830 | 4592 | 3411 | 1.3 |
View observed and expected dry season malaria cases
# Helper function to create maps of observed and expected dry season malaria cases
create.malaria.map <- function(malaria.data,
variable = NA,
title = NA,
legend.title = NA){
# observed and expected malaria incidence map
# malaria.data: data frame containing observed and expected malaria cases
# variable: variable name (as character, in quotes e.g. "observed")
# title: map title in quotes
# legend.title: legend title in quotes
# returns:
# a tmap-element (plots a map)
tm_shape(malaria.data)+
tm_fill(col = variable,
breaks = c(0, 500, 1000, 2500, 5000, 10000, 15000, 20000, 25000),
palette = "YlOrRd",
title = legend.title)+
tm_borders(lw = 0.3)+
tm_layout(legend.position = c(0.75,"bottom"),
legend.text.size = 0.5,
legend.title.size = 0.7,
frame = FALSE)+
tm_credits(title,
position = c(0.2, 0.8),
size = 1)
}
# Invoking the function
# 2017 observed and expected malaria cases -------------------------------------
observed_malaria_2017_map <- create.malaria.map(malaria_pop_by_catchment_2017,
variable = "dr_2017",
title = "2017",
legend.title = "Observed malaria")
expected_malaria_2017_map <- create.malaria.map(expected_malaria_2017,
variable = "expected_2017",
title = "2017",
legend.title = "Expected malaria")
# 2018 observed and expected malaria cases -------------------------------------
observed_malaria_2018_map <- create.malaria.map(malaria_pop_by_catchment_2018,
variable = "dr_2018",
title = "2018",
legend.title = "Observed malaria")
expected_malaria_2018_map <- create.malaria.map(expected_malaria_2018,
variable = "expected_2018",
title = "2018",
legend.title = "Expected malaria")
# 2019 observed and expected malaria cases -------------------------------------
observed_malaria_2019_map <- create.malaria.map(malaria_pop_by_catchment_2019,
variable = "dr_2019",
title = "2019",
legend.title = "Observed malaria")
expected_malaria_2019_map <- create.malaria.map(expected_malaria_2019,
variable = "expected_2019",
title = "2019",
legend.title = "Expected malaria")
# 2020 observed and expected malaria cases -------------------------------------
observed_malaria_2020_map <- create.malaria.map(malaria_pop_by_catchment_2020,
variable = "dr_2020",
title = "2020",
legend.title = "Observed malaria")
expected_malaria_2020_map <- create.malaria.map(expected_malaria_2020,
variable = "expected_2020",
title = "2020",
legend.title = "Expected malaria")
# Layout maps ------------------------------------------------------------------
tmap::tmap_arrange(observed_malaria_2017_map, expected_malaria_2017_map,
observed_malaria_2018_map, expected_malaria_2018_map,
observed_malaria_2019_map, expected_malaria_2019_map,
observed_malaria_2020_map, expected_malaria_2020_map, ncol = 2)
Fig 6: Observed and expected malaria incidence by health facility catchment area, Kasungu
View SMR by catchment
A ratio greater than 1.0 indicates that more malaria cases have occurred than would have been expected, while a ratio less than 1.0 indicates that less cases have occurred.
# max(SMR_2017$SMR)
# [1] 6.93
# > max(SMR_2018$SMR)
# [1] 7.93
# > max(SMR_2019$SMR)
# [1] 7.71
# > max(SMR_2020$SMR)
# [1] 5.68
# Helper function to create maps of SMR by catchment ----------------------------------
create.smr.map <- function(smr.data,
variable = "SMR",
title = NA,
legend.title = "SMR"){
# SMR by catchment map
# smr.data: sf polygon object containing SMR by catchment data
# variable: variable name (as character, in qoutes)
# title: map title in quotes
# legend.title: legend title in qoutes
# returns:
# a tmap-element (plots a map)
tm_shape(smr.data)+
tm_fill(col = variable,
breaks = c(0, 0.5, 1, 1.5, 2, 2.5, 5, 8),
palette = "-magma",
title = legend.title)+
tm_borders(lw = 0.3)+
tm_layout(legend.position = c(0.75,"bottom"),
legend.text.size = 0.5,
legend.title.size = 0.7,
frame = FALSE)+
tm_credits(title,
position = c(0.2, 0.8),
size = 1)
}
# Invoking function -------------------------------------------------------------------
SMR_2017_map <- create.smr.map(SMR_2017, title = "2017")
SMR_2018_map <- create.smr.map(SMR_2018, title = "2018")
SMR_2019_map <- create.smr.map(SMR_2019, title = "2019")
SMR_2020_map <- create.smr.map(SMR_2020, title = "2020")
# Layout maps -------------------------------------------------------------------------
tmap::tmap_arrange(SMR_2017_map, SMR_2018_map, SMR_2019_map, SMR_2020_map, ncol = 2)
Fig. 7: Standardised morbidity ratio of malaria by health facility catchment
Calculate the proportion of the catchment population living within 1km, 2km, 3km of water bodies
First, using st_buffer, we compute 1km, 2km and 3km buffers around dry season water bodies obtained from LandSat satellite imagery using TropWet tool in Google Earth Engine. Then geometry of the buffer features are then combined resulting in resolved internal boundaries to enable extracting population values from WorldPop raster. Finally, we calculate the proportion of people in each catchment area living within water bodies.
# Combine and transform TropWet derived waterbody polygons -------------------------------
surface_waterbodies_2017 <- sf::st_as_sf(
st_cast(
st_union(
st_buffer(dryseason_waterbodies_2017, dist = 30)), "POLYGON"))
surface_waterbodies_2018 <- sf::st_as_sf(
st_cast(
st_union(
st_buffer(dryseason_waterbodies_2018, dist = 30)), "POLYGON"))
surface_waterbodies_2019 <- sf::st_as_sf(
st_cast(
st_union(
st_buffer(dryseason_waterbodies_2019, dist = 30)), "POLYGON"))
surface_waterbodies_2020 <- sf::st_as_sf(
st_cast(
st_union(
st_buffer(dryseason_waterbodies_2020, dist = 30)), "POLYGON"))
# Helper function to compute 1km, 2km and 3km buffers around the water bodies ---------------------
create.waterbody.buffer <- function(waterbody, distance, catchment){
# function for creating buffers around waterbodies
# arguments:
# waterbody: waterbody shapefile
# distance: buffer distance in meters
# catchment: catchment area shapefile
# returns:
# buffered waterbodies
# Create buffers around water bodies
buffer_radius <- sf::st_buffer(waterbody, distance)
# Dissolve the buffers
buffer_union <- sf::st_as_sf(st_cast(st_union(buffer_radius),"MULTIPOLYGON"))
# Assign attributes of the 'catchment' to each of the water bodies.
buffer_intersect <- sf::st_intersection(buffer_union, catchment)
buffer_intersect_sf <- sf::st_as_sf(buffer_intersect)
# Convert the MULTIPOLYGON object into several POLYGON objects
buffer_intersect_polygons <- sf::st_cast(
sf::st_buffer(buffer_intersect_sf,0.0), "MULTIPOLYGON") %>%
sf::st_cast("POLYGON")
# Polygons being seen to be in multiple catchments
sf::st_intersects(buffer_intersect_polygons, catchment)
# Make the assumption that the attribute is constant throughout the geometry
sf::st_agr(buffer_intersect_polygons) = "constant"
sf::st_agr(catchment) = "constant"
return(out = buffer_intersect_polygons)
}
# Invoking function
# For 2017 TropWet surface water polygons --------------------------------------------------------
buffer_1km_2017 <- create.waterbody.buffer(waterbody = surface_waterbodies_2017,
distance = 1000,
catchment = malire_new)
buffer_2km_2017 <- create.waterbody.buffer(waterbody = surface_waterbodies_2017,
distance = 2000,
catchment = malire_new)
buffer_3km_2017 <- create.waterbody.buffer(waterbody = surface_waterbodies_2017,
distance = 3000,
catchment = malire_new)
# For 2018 TropWet surface water polygons --------------------------------------------------------
buffer_1km_2018 <- create.waterbody.buffer(waterbody = surface_waterbodies_2018,
distance = 1000,
catchment = malire_new)
buffer_2km_2018 <- create.waterbody.buffer(waterbody = surface_waterbodies_2018,
distance = 2000,
catchment = malire_new)
buffer_3km_2018 <- create.waterbody.buffer(waterbody = surface_waterbodies_2018,
distance = 3000,
catchment = malire_new)
# For 2019 TropWet surface water polygons ------------------------------------------------------
buffer_1km_2019 <- create.waterbody.buffer(waterbody = surface_waterbodies_2019,
distance = 1000,
catchment = malire_new)
buffer_2km_2019 <- create.waterbody.buffer(waterbody = surface_waterbodies_2019,
distance = 2000,
catchment = malire_new)
buffer_3km_2019 <- create.waterbody.buffer(waterbody = surface_waterbodies_2019,
distance = 3000,
catchment = malire_new)
# For 2020 TropWet surface water polygons ------------------------------------------------------
buffer_1km_2020 <- create.waterbody.buffer(waterbody = surface_waterbodies_2020,
distance = 1000,
catchment = malire_new)
buffer_2km_2020 <- create.waterbody.buffer(waterbody = surface_waterbodies_2020,
distance = 2000,
catchment = malire_new)
buffer_3km_2020 <- create.waterbody.buffer(waterbody = surface_waterbodies_2020,
distance = 3000,
catchment = malire_new)View the created waterbody buffers
# Map the buffers
create.buffer.map <- function(buffers, boundary = malire_new, title = NA){
# function for creating buffer map in ggplot
# arguments:
# buffer: waterbodies buffer polygon layer
# boundary: health facility catchment polygons
# title: main title
# returns:
# a map-element (plots a map)
ggplot(data = buffers)+
geom_sf()+
geom_sf(data = boundary,
fill = NA)+
theme_void()+
labs(title = title)
}
# Invoking the function
# For 2017 -------------------------------------------------------------------------------
buffer_1km_2017_map <- create.buffer.map(buffer_1km_2017, title = "2017: 1km Buffers")
buffer_2km_2017_map <- create.buffer.map(buffer_2km_2017, title = "2017: 2km Buffers")
buffer_3km_2017_map <- create.buffer.map(buffer_3km_2017, title = "2017: 3km Buffers")
# For 2018 --------------------------------------------------------------------------------
buffer_1km_2018_map <- create.buffer.map(buffer_1km_2018, title = "2018: 1km Buffers")
buffer_2km_2018_map <- create.buffer.map(buffer_2km_2018, title = "2018: 2km Buffers")
buffer_3km_2018_map <- create.buffer.map(buffer_3km_2018, title = "2018: 3km Buffers")
# For 2019 ---------------------------------------------------------------------------------
buffer_1km_2019_map <- create.buffer.map(buffer_1km_2019, title = "2019: 1km Buffers")
buffer_2km_2019_map <- create.buffer.map(buffer_2km_2019, title = "2019: 2km Buffers")
buffer_3km_2019_map <- create.buffer.map(buffer_3km_2019, title = "2019: 3km Buffers")
# For 2020 --------------------------------------------------------------------------------
buffer_1km_2020_map <- create.buffer.map(buffer_1km_2020, title = "2020: 1km Buffers")
buffer_2km_2020_map <- create.buffer.map(buffer_2km_2020, title = "2020: 2km Buffers")
buffer_3km_2020_map <- create.buffer.map(buffer_3km_2020, title = "2020: 3km Buffers")
grid.arrange(buffer_1km_2017_map, buffer_1km_2018_map, buffer_1km_2019_map, buffer_1km_2020_map,
buffer_2km_2017_map, buffer_2km_2018_map, buffer_2km_2019_map, buffer_2km_2020_map,
buffer_3km_2017_map, buffer_3km_2018_map, buffer_3km_2019_map, buffer_3km_2020_map, ncol = 4)
Fig 8. Buffers around dry season waterbodies in Kasungu
Extract the population living within waterbody buffers by catchment area
# Helper function to calculate estimated number of people living within waterbody buffers
# in each catchment area
estimate.buffer.pop <- function(catchment.population, buffers, catchment.area){
# Extract population estimates from WorldPop raster
buffers$buffer_pop <- raster::extract(catchment.population,
buffers,
fun = sum,
na.rm = TRUE)
# Find which catchment each polygon belongs to using its centroid - a point dataset
# representing the geographic center-points of the polygons
buffer_by_catchment <- st_intersection(st_centroid(buffers), catchment.area)
# Notice that the buffer_catchment is comprised of separate POLYGONS (buffer_by_catchment$x).
# The first step is to “dissolve” away these POLYGONS into one MULTIPOLYGON.
# There is no sf equivalent to the QGIS or ArcMap “dissolve” operation.
# Instead we use a combination of group_by and summarize from the dplyr package.
# Stats::aggregate from sf package, and dplyr::summarize both do essentially the same.
buffer_pop_aggregated <- buffer_by_catchment %>%
dplyr::group_by(DN) %>%
dplyr::summarize(
buffer_pop_aggregated = round(sum(buffer_pop, na.rm = TRUE)))
buffer_pop <- merge(
catchment.area, st_drop_geometry(
buffer_pop_aggregated), by = 'DN', all.x = TRUE)
return(out = buffer_pop)
}
# Invoking the function and calculating proportion of
# catchment population living within buffers
# 2017 buffer population -------------------------------------------------------
buffer_pop_1km_2017 <- estimate.buffer.pop(
kasungu_population_2017,
buffer_1km_2017,
malaria_pop_by_catchment_2017) %>%
dplyr::rename(catchment_pop = pop,
buffer_pop = buffer_pop_aggregated) %>%
dplyr::mutate(
prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))
buffer_pop_2km_2017 <- estimate.buffer.pop(
kasungu_population_2017,
buffer_2km_2017,
malaria_pop_by_catchment_2017) %>%
dplyr::rename(catchment_pop = pop,
buffer_pop = buffer_pop_aggregated) %>%
dplyr::mutate(
prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))
buffer_pop_3km_2017 <- estimate.buffer.pop(
kasungu_population_2017,
buffer_3km_2017,
malaria_pop_by_catchment_2017) %>%
dplyr::rename(catchment_pop = pop,
buffer_pop = buffer_pop_aggregated) %>%
dplyr::mutate(
prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))
# 2018 buffer population -------------------------------------------------------
buffer_pop_1km_2018 <- estimate.buffer.pop(
kasungu_population_2018,
buffer_1km_2018,
malaria_pop_by_catchment_2018) %>%
dplyr::rename(catchment_pop = pop,
buffer_pop = buffer_pop_aggregated) %>%
dplyr::mutate(
prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))
buffer_pop_2km_2018 <- estimate.buffer.pop(
kasungu_population_2018,
buffer_2km_2018,
malaria_pop_by_catchment_2018) %>%
dplyr::rename(catchment_pop = pop,
buffer_pop = buffer_pop_aggregated) %>%
dplyr::mutate(
prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))
buffer_pop_3km_2018 <- estimate.buffer.pop(
kasungu_population_2018,
buffer_3km_2018,
malaria_pop_by_catchment_2018) %>%
dplyr::rename(catchment_pop = pop,
buffer_pop = buffer_pop_aggregated) %>%
dplyr::mutate(
prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))
# 2019 buffer population -------------------------------------------------------
buffer_pop_1km_2019 <- estimate.buffer.pop(
kasungu_population_2019,
buffer_1km_2019,
malaria_pop_by_catchment_2019) %>%
dplyr::rename(catchment_pop = pop,
buffer_pop = buffer_pop_aggregated) %>%
dplyr::mutate(
prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))
buffer_pop_2km_2019 <- estimate.buffer.pop(
kasungu_population_2019,
buffer_2km_2019,
malaria_pop_by_catchment_2019) %>%
dplyr::rename(catchment_pop = pop,
buffer_pop = buffer_pop_aggregated) %>%
dplyr::mutate(
prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))
buffer_pop_3km_2019 <- estimate.buffer.pop(
kasungu_population_2019,
buffer_3km_2019,
malaria_pop_by_catchment_2019) %>%
dplyr::rename(catchment_pop = pop,
buffer_pop = buffer_pop_aggregated) %>%
dplyr::mutate(
prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))
# 2020 buffer population -------------------------------------------------------
buffer_pop_1km_2020 <- estimate.buffer.pop(
kasungu_population_2020,
buffer_1km_2020,
malaria_pop_by_catchment_2020) %>%
dplyr::rename(catchment_pop = pop,
buffer_pop = buffer_pop_aggregated) %>%
dplyr::mutate(
prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))
buffer_pop_2km_2020 <- estimate.buffer.pop(
kasungu_population_2020,
buffer_2km_2020,
malaria_pop_by_catchment_2020) %>%
dplyr::rename(catchment_pop = pop,
buffer_pop = buffer_pop_aggregated) %>%
dplyr::mutate(
prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))
buffer_pop_3km_2020 <- estimate.buffer.pop(
kasungu_population_2020,
buffer_3km_2020,
malaria_pop_by_catchment_2020) %>%
dplyr::rename(catchment_pop = pop,
buffer_pop = buffer_pop_aggregated) %>%
dplyr::mutate(
prop_buffer_catchment_pop = round((buffer_pop/catchment_pop)*100))Mapping proportion of catchment population living within waterbodies
# Helper function to create maps of proportion of people living in proximity ----------
# to water bodies in each catchment area
create.pop.proportion.map <- function(pop.data,
variable = "prop_buffer_catchment_pop",
title = NA,
legend.title = NA){
# pop.data: sf polygon object containing proportion of catchment population
# living within water bodies
# variable: variable name (as character, in qoutes)
# title: map title in quotes
# legend.title: legend title in qoutes
# returns:
# a tmap-element (plots a map)
tm_shape(pop.data)+
tm_fill(col = variable,
breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100),
palette = "YlOrBr",
title = legend.title)+
tm_borders(lw = 0.3)+
tm_layout(legend.position = c(0.8,"bottom"),
legend.text.size = 0.5,
legend.title.size = 0.7,
frame = FALSE)+
tm_credits(title,
position = c(0.25, 0.75),
size = 1)
}
# Invoking function
# 2017 population proportion ---------------------------------------------------
pop_proportion_1km_2017_map <- create.pop.proportion.map(
buffer_pop_1km_2017,
title = "2017",
legend.title = "Population within \n1km buffers (%)")
pop_proportion_2km_2017_map <- create.pop.proportion.map(
buffer_pop_2km_2017,
title = "2017",
legend.title = "Population within \n2km buffers (%)")
pop_proportion_3km_2017_map <- create.pop.proportion.map(
buffer_pop_3km_2017,
title = "2017",
legend.title = "Population within \n3km buffers (%)")
# 2018 population proportion ---------------------------------------------------
pop_proportion_1km_2018_map <- create.pop.proportion.map(
buffer_pop_1km_2018,
title = "2018",
legend.title = "Population within \n1km buffers (%)")
pop_proportion_2km_2018_map <- create.pop.proportion.map(
buffer_pop_2km_2018,
title = "2018",
legend.title = "Population within \n2km buffers (%)")
pop_proportion_3km_2018_map <- create.pop.proportion.map(
buffer_pop_3km_2018,
title = "2018",
legend.title = "Population within \n3km buffers (%)")
# 2019 population proportion ---------------------------------------------------
pop_proportion_1km_2019_map <- create.pop.proportion.map(
buffer_pop_1km_2019,
title = "2019",
legend.title = "Population within \n1km buffers (%)")
pop_proportion_2km_2019_map <- create.pop.proportion.map(
buffer_pop_2km_2019,
title = "2019",
legend.title = "Population within \n2km buffers (%)")
pop_proportion_3km_2019_map <- create.pop.proportion.map(
buffer_pop_3km_2019,
title = "2019",
legend.title = "Population within \n3km buffers (%)")
# 2020 population proportion ---------------------------------------------------
pop_proportion_1km_2020_map <- create.pop.proportion.map(
buffer_pop_1km_2020,
title = "2020",
legend.title = "Population within \n1km buffers (%)")
pop_proportion_2km_2020_map <- create.pop.proportion.map(
buffer_pop_2km_2020,
title = "2020",
legend.title = "Population within \n2km buffers (%)")
pop_proportion_3km_2020_map <- create.pop.proportion.map(
buffer_pop_3km_2020,
title = "2020",
legend.title = "Population within \n3km buffers (%)")
# Layout maps ------------------------------------------------------------------
tmap::tmap_arrange(pop_proportion_1km_2017_map, pop_proportion_2km_2017_map,
pop_proportion_3km_2017_map, pop_proportion_1km_2018_map,
pop_proportion_2km_2018_map, pop_proportion_3km_2018_map,
pop_proportion_1km_2019_map, pop_proportion_2km_2019_map,
pop_proportion_3km_2019_map, pop_proportion_1km_2020_map,
pop_proportion_2km_2020_map, pop_proportion_3km_2020_map, ncol = 3)
Fig 8. Proportion of catchment population living around water bodies
Scatter plots of SMR against the proportion of the catchment population living waterbody buffers
A correlation coeeficient of more than zero (cor.coeff r > 0.1) indicates some positive association between the SMR and the buffer population variables. That is, SMR of dry season malaria increases with increase in number of people suurounding water bodies.
# Helper function to tidy and bind the SMR and proportion of -------------------
# buffer-catchment population data frames
tidy.data <- function(smr.df,
proportion.pop.1km,
proprotion.pop.2km,
proportion.pop.3km){
# Convert the sf object to dataframes-------------------------------------------
smr_df <- as.data.frame(smr.df) %>%
dplyr::select(rowID, Names, SMR)
proportion_pop_1km_df <- as.data.frame(proportion.pop.1km) %>%
dplyr::select(rowID, prop_pop_1km = `prop_buffer_catchment_pop`)
proportion_pop_2km_df <- as.data.frame(proprotion.pop.2km)%>%
dplyr::select(rowID, prop_pop_2km = `prop_buffer_catchment_pop`)
proportion_pop_3km_df <- as.data.frame(proportion.pop.3km)%>%
dplyr::select(rowID, prop_pop_3km = `prop_buffer_catchment_pop`)
# Merge SMR and population data frames -----------------------------------------
combined_1 <- merge(smr_df, proportion_pop_1km_df, by = "rowID", all = TRUE)
combined_2 <- merge(proportion_pop_2km_df, proportion_pop_3km_df)
combined_fully <- merge(combined_1, combined_2, by = "rowID", all = TRUE)
}
# Invoking the function --------------------------------------------------------
smr_pop_2017 <- tidy.data(SMR_2017, buffer_pop_1km_2017, buffer_pop_2km_2017, buffer_pop_3km_2017)
smr_pop_2018 <- tidy.data(SMR_2018, buffer_pop_1km_2018, buffer_pop_2km_2018, buffer_pop_3km_2018)
smr_pop_2019 <- tidy.data(SMR_2019, buffer_pop_1km_2019, buffer_pop_2km_2019, buffer_pop_3km_2019)
smr_pop_2020 <- tidy.data(SMR_2020, buffer_pop_1km_2020, buffer_pop_2km_2020, buffer_pop_3km_2020)
# Helper function to create scatter plots --------------------------------------
create.scatter.plot <- function(smr.pop.df,
independent.var = NA,
dependent.var = "SMR",
x.label = NA,
plot.title = NA){
scatter.plot <- ggpubr::ggscatter(smr.pop.df, # data frame
x = independent.var, # x-axis variable
y = dependent.var, # y-axis variable
add = "reg.line", # Add regression line
conf.int = TRUE, # Add confidence interval
add.params = list(color = "red",
fill = "lightgray"),
palette = "jco", # journal color palette. see ?ggpar
xlab = x.label, # x-axis label
ylab = "SMR", # y-axis label
title = plot.title)+
ggpubr::stat_cor(label.y = 5)+ # Add correlation coefficient
ggpubr::font("title", size = 10, face = "bold")+
ggpubr::font("xlab", size = 10)+
ggpubr::font("ylab", size = 10)
return(scatter.plot)
}
# Invoking function
# 2017 scatter plots ------------------------------------------------------------
scatter_1km_2017 <- create.scatter.plot(smr_pop_2017, independent.var = "prop_pop_1km",
x.label = "Percentage of catchment population \nliving in 1km buffer",
plot.title = "2017")
scatter_2km_2017 <- create.scatter.plot(smr_pop_2017, independent.var = "prop_pop_2km",
x.label = "Percentage of catchment population \nliving in 2km buffer",
plot.title = "2017")
scatter_3km_2017 <- create.scatter.plot(smr_pop_2017, independent.var = "prop_pop_3km",
x.label = "Percentage of catchment population \nliving in 3km buffer",
plot.title = "2017")
# 2018 scatter plots -----------------------------------------------------------
scatter_1km_2018 <- create.scatter.plot(smr_pop_2018, independent.var = "prop_pop_1km",
x.label = "Percentage of catchment population \nliving in 1km buffer",
plot.title = "2018")
scatter_2km_2018 <- create.scatter.plot(smr_pop_2018, independent.var = "prop_pop_2km",
x.label = "Percentage of catchment population \nliving in 2km buffer",
plot.title = "2018")
scatter_3km_2018 <- create.scatter.plot(smr_pop_2018, independent.var = "prop_pop_3km",
x.label = "Percentage of catchment population \nliving in 3km buffer",
plot.title = "2018")
# 2019 scatter plots -----------------------------------------------------------
scatter_1km_2019 <- create.scatter.plot(smr_pop_2019, independent.var = "prop_pop_1km",
x.label = "Percentage of catchment population \nliving in 1km buffer",
plot.title = "2019")
scatter_2km_2019 <- create.scatter.plot(smr_pop_2019, independent.var = "prop_pop_2km",
x.label = "Percentage of catchment population \nliving in 2km buffer",
plot.title = "2019")
scatter_3km_2019 <- create.scatter.plot(smr_pop_2019, independent.var = "prop_pop_3km",
x.label = "Percentage of catchment population \nliving in 3km buffer",
plot.title = "2019")
# 2020 scatter plots -----------------------------------------------------------
scatter_1km_2020 <- create.scatter.plot(smr_pop_2020, independent.var = "prop_pop_1km",
x.label = "Percentage of catchment population \nliving in 1km buffer",
plot.title = "2020")
scatter_2km_2020 <- create.scatter.plot(smr_pop_2020, independent.var = "prop_pop_2km",
x.label = "Percentage of catchment population \nliving in 2km buffer",
plot.title = "2020")
scatter_3km_2020 <- create.scatter.plot(smr_pop_2020, independent.var = "prop_pop_3km",
x.label = "Percentage of catchment population \nliving in 3km buffer",
plot.title = "2020")
# Arrange the plots ------------------------------------------------------------
ggpubr::ggarrange(scatter_1km_2017, scatter_2km_2017, scatter_3km_2017,
scatter_1km_2018, scatter_2km_2018, scatter_3km_2019,
scatter_1km_2019, scatter_2km_2019, scatter_3km_2019,
scatter_1km_2020, scatter_2km_2020, scatter_3km_2020,
ncol = 3, nrow = 4)
Fig 10. Relationship between standardised morbidity ratio and living near waterbodies
Model fitting
# Combine data for model fitting -----------------------------------------------
model_data_2017 <- merge(expected_malaria_2017, smr_pop_2017, by = "rowID", all = TRUE) %>%
dplyr::select(-Names.y) %>%
dplyr::rename(Names = Names.x)
model_data_2018 <- merge(expected_malaria_2018, smr_pop_2018, by = "rowID", all = TRUE) %>%
dplyr::select(-Names.y) %>%
dplyr::rename(Names = Names.x)
model_data_2019 <- merge(expected_malaria_2019, smr_pop_2019, by = "rowID", all = TRUE) %>%
dplyr::select(-Names.y) %>%
dplyr::rename(Names = Names.x)
model_data_2020 <- merge(expected_malaria_2020, smr_pop_2020, by = "rowID", all = TRUE) %>%
dplyr::select(-Names.y) %>%
dplyr::rename(Names = Names.x)
# Fit generalised linear model -------------------------------------------------
# Defining model parameters:
# response variable: observed_2017, observed_2018, observed_2019, observed_2020 are
# recorded dry season malaria cases in that year
# risk factor: prop_pop_1km, prop_pop_2km, prop_pop_3km are the percentage of people living
# within 1km, 2km and 3km buffers of water bodies, respectively.
# offset: expected_year is the number of malaria cases we would expect if the malaria rate
# was equal in all the catchment areas
# 2017 -------------------------------------------------------------------------
model_1km_2017 <- glm(observed_2017~1+prop_pop_1km+offset(log(expected_2017)),
data = model_data_2017, family = 'poisson')
summary(model_1km_2017)##
## Call:
## glm(formula = observed_2017 ~ 1 + prop_pop_1km + offset(log(expected_2017)),
## family = "poisson", data = model_data_2017)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -50.150 -19.463 -1.324 20.443 118.758
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.2223780 0.0059167 -37.58 <2e-16 ***
## prop_pop_1km 0.0247164 0.0004965 49.78 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 33562 on 23 degrees of freedom
## Residual deviance: 31130 on 22 degrees of freedom
## (3 observations deleted due to missingness)
## AIC: 31372
##
## Number of Fisher Scoring iterations: 5sjPlot::tab_model(model_1km_2017)| observed_2017 | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| (Intercept) | 0.80 | 0.79 – 0.81 | <0.001 |
| prop_pop_1km | 1.03 | 1.02 – 1.03 | <0.001 |
| Observations | 24 | ||
| R2 Nagelkerke | 1.000 | ||
model_2km_2017 <- glm(observed_2017~1+prop_pop_2km+offset(log(expected_2017)),
data = model_data_2017, family = 'poisson')
summary(model_2km_2017)##
## Call:
## glm(formula = observed_2017 ~ 1 + prop_pop_2km + offset(log(expected_2017)),
## family = "poisson", data = model_data_2017)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -45.875 -19.972 1.747 19.501 118.324
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3126953 0.0062855 -49.75 <2e-16 ***
## prop_pop_2km 0.0134272 0.0002117 63.42 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 33564 on 24 degrees of freedom
## Residual deviance: 29695 on 23 degrees of freedom
## (2 observations deleted due to missingness)
## AIC: 29947
##
## Number of Fisher Scoring iterations: 5sjPlot::tab_model(model_2km_2017)| observed_2017 | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| (Intercept) | 0.73 | 0.72 – 0.74 | <0.001 |
| prop_pop_2km | 1.01 | 1.01 – 1.01 | <0.001 |
| Observations | 25 | ||
| R2 Nagelkerke | 1.000 | ||
model_3km_2017 <- glm(observed_2017~1+prop_pop_3km+offset(log(expected_2017)),
data = model_data_2017, family = 'poisson')
summary(model_3km_2017)##
## Call:
## glm(formula = observed_2017 ~ 1 + prop_pop_3km + offset(log(expected_2017)),
## family = "poisson", data = model_data_2017)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -47.915 -21.148 1.191 19.098 121.065
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3936999 0.0075710 -52.00 <2e-16 ***
## prop_pop_3km 0.0100819 0.0001629 61.88 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 33564 on 24 degrees of freedom
## Residual deviance: 29726 on 23 degrees of freedom
## (2 observations deleted due to missingness)
## AIC: 29977
##
## Number of Fisher Scoring iterations: 5sjPlot::tab_model(model_3km_2017)| observed_2017 | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| (Intercept) | 0.67 | 0.66 – 0.68 | <0.001 |
| prop_pop_3km | 1.01 | 1.01 – 1.01 | <0.001 |
| Observations | 25 | ||
| R2 Nagelkerke | 1.000 | ||
# 2018 -------------------------------------------------------------------------
model_1km_2018 <- glm(observed_2018~1+prop_pop_1km+offset(log(expected_2018)),
data = model_data_2018, family = 'poisson')
summary(model_1km_2018)##
## Call:
## glm(formula = observed_2018 ~ 1 + prop_pop_1km + offset(log(expected_2018)),
## family = "poisson", data = model_data_2018)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -55.722 -12.556 3.871 13.940 137.674
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3987643 0.0060731 -65.66 <2e-16 ***
## prop_pop_1km 0.0323262 0.0004041 79.99 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 42339 on 25 degrees of freedom
## Residual deviance: 35921 on 24 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 36187
##
## Number of Fisher Scoring iterations: 5sjPlot::tab_model(model_1km_2018)| observed_2018 | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| (Intercept) | 0.67 | 0.66 – 0.68 | <0.001 |
| prop_pop_1km | 1.03 | 1.03 – 1.03 | <0.001 |
| Observations | 26 | ||
| R2 Nagelkerke | 1.000 | ||
model_2km_2018 <- glm(observed_2018~1+prop_pop_2km+offset(log(expected_2018)),
data = model_data_2018, family = 'poisson')
summary(model_2km_2018)##
## Call:
## glm(formula = observed_2018 ~ 1 + prop_pop_2km + offset(log(expected_2018)),
## family = "poisson", data = model_data_2018)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -58.860 -11.691 4.018 20.409 114.785
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5633136 0.0067344 -83.65 <2e-16 ***
## prop_pop_2km 0.0167495 0.0001679 99.76 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 42339 on 25 degrees of freedom
## Residual deviance: 32218 on 24 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 32484
##
## Number of Fisher Scoring iterations: 5sjPlot::tab_model(model_2km_2018)| observed_2018 | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| (Intercept) | 0.57 | 0.56 – 0.58 | <0.001 |
| prop_pop_2km | 1.02 | 1.02 – 1.02 | <0.001 |
| Observations | 26 | ||
| R2 Nagelkerke | 1.000 | ||
model_3km_2018 <- glm(observed_2018~1+prop_pop_3km+offset(log(expected_2018)),
data = model_data_2018, family = 'poisson')
summary(model_3km_2018)##
## Call:
## glm(formula = observed_2018 ~ 1 + prop_pop_3km + offset(log(expected_2018)),
## family = "poisson", data = model_data_2018)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -52.920 -14.159 5.671 22.270 106.698
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8634330 0.0089250 -96.74 <2e-16 ***
## prop_pop_3km 0.0165420 0.0001522 108.66 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 42339 on 25 degrees of freedom
## Residual deviance: 29688 on 24 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 29955
##
## Number of Fisher Scoring iterations: 5sjPlot::tab_model(model_3km_2018)| observed_2018 | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| (Intercept) | 0.42 | 0.41 – 0.43 | <0.001 |
| prop_pop_3km | 1.02 | 1.02 – 1.02 | <0.001 |
| Observations | 26 | ||
| R2 Nagelkerke | 1.000 | ||
# 2019 -------------------------------------------------------------------------
model_1km_2019 <- glm(observed_2019~1+prop_pop_1km+offset(log(expected_2019)),
data = model_data_2019, family = 'poisson')
summary(model_1km_2019)##
## Call:
## glm(formula = observed_2019 ~ 1 + prop_pop_1km + offset(log(expected_2019)),
## family = "poisson", data = model_data_2019)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -66.815 -18.307 1.353 23.118 102.556
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3353561 0.0077298 -43.38 <2e-16 ***
## prop_pop_1km 0.0180337 0.0003618 49.84 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 33446 on 26 degrees of freedom
## Residual deviance: 30926 on 25 degrees of freedom
## AIC: 31194
##
## Number of Fisher Scoring iterations: 5sjPlot::tab_model(model_1km_2019)| observed_2019 | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| (Intercept) | 0.72 | 0.70 – 0.73 | <0.001 |
| prop_pop_1km | 1.02 | 1.02 – 1.02 | <0.001 |
| Observations | 27 | ||
| R2 Nagelkerke | 1.000 | ||
model_2km_2019 <- glm(observed_2019~1+prop_pop_2km+offset(log(expected_2019)),
data = model_data_2019, family = 'poisson')
summary(model_2km_2019)##
## Call:
## glm(formula = observed_2019 ~ 1 + prop_pop_2km + offset(log(expected_2019)),
## family = "poisson", data = model_data_2019)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -62.833 -21.458 2.378 24.375 101.995
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3051901 0.0088349 -34.54 <2e-16 ***
## prop_pop_2km 0.0066732 0.0001748 38.17 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 33446 on 26 degrees of freedom
## Residual deviance: 31955 on 25 degrees of freedom
## AIC: 32223
##
## Number of Fisher Scoring iterations: 5sjPlot::tab_model(model_2km_2019)| observed_2019 | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| (Intercept) | 0.74 | 0.72 – 0.75 | <0.001 |
| prop_pop_2km | 1.01 | 1.01 – 1.01 | <0.001 |
| Observations | 27 | ||
| R2 Nagelkerke | 1.000 | ||
model_3km_2019 <- glm(observed_2019~1+prop_pop_3km+offset(log(expected_2019)),
data = model_data_2019, family = 'poisson')
summary(model_3km_2019)##
## Call:
## glm(formula = observed_2019 ~ 1 + prop_pop_3km + offset(log(expected_2019)),
## family = "poisson", data = model_data_2019)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -60.344 -22.596 4.515 23.613 102.074
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3214559 0.0105946 -30.34 <2e-16 ***
## prop_pop_3km 0.0046685 0.0001437 32.49 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 33446 on 26 degrees of freedom
## Residual deviance: 32357 on 25 degrees of freedom
## AIC: 32624
##
## Number of Fisher Scoring iterations: 5sjPlot::tab_model(model_3km_2019)| observed_2019 | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| (Intercept) | 0.73 | 0.71 – 0.74 | <0.001 |
| prop_pop_3km | 1.00 | 1.00 – 1.00 | <0.001 |
| Observations | 27 | ||
| R2 Nagelkerke | 1.000 | ||
# 2020 -------------------------------------------------------------------------
model_1km_2020 <- glm(observed_2020~1+prop_pop_1km+offset(log(expected_2020)),
data = model_data_2020, family = 'poisson')
summary(model_1km_2020)##
## Call:
## glm(formula = observed_2020 ~ 1 + prop_pop_1km + offset(log(expected_2020)),
## family = "poisson", data = model_data_2020)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -83.14 -23.49 9.60 24.68 89.43
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.0236328 0.0060477 -3.908 9.32e-05 ***
## prop_pop_1km 0.0024199 0.0004792 5.050 4.41e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 36512 on 22 degrees of freedom
## Residual deviance: 36486 on 21 degrees of freedom
## (4 observations deleted due to missingness)
## AIC: 36719
##
## Number of Fisher Scoring iterations: 5sjPlot::tab_model(model_1km_2020)| observed_2020 | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| (Intercept) | 0.98 | 0.97 – 0.99 | <0.001 |
| prop_pop_1km | 1.00 | 1.00 – 1.00 | <0.001 |
| Observations | 23 | ||
| R2 Nagelkerke | 0.670 | ||
model_2km_2020 <- glm(observed_2020~1+prop_pop_2km+offset(log(expected_2020)),
data = model_data_2020, family = 'poisson')
summary(model_2km_2020)##
## Call:
## glm(formula = observed_2020 ~ 1 + prop_pop_2km + offset(log(expected_2020)),
## family = "poisson", data = model_data_2020)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -79.646 -30.318 7.777 26.260 87.868
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.1690115 0.0065167 -25.93 <2e-16 ***
## prop_pop_2km 0.0056305 0.0002104 26.76 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 37822 on 23 degrees of freedom
## Residual deviance: 37108 on 22 degrees of freedom
## (3 observations deleted due to missingness)
## AIC: 37350
##
## Number of Fisher Scoring iterations: 5sjPlot::tab_model(model_2km_2020)| observed_2020 | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| (Intercept) | 0.84 | 0.83 – 0.86 | <0.001 |
| prop_pop_2km | 1.01 | 1.01 – 1.01 | <0.001 |
| Observations | 24 | ||
| R2 Nagelkerke | 1.000 | ||
model_3km_2020 <- glm(observed_2020~1+prop_pop_3km+offset(log(expected_2020)),
data = model_data_2020, family = 'poisson')
summary(model_3km_2020)##
## Call:
## glm(formula = observed_2020 ~ 1 + prop_pop_3km + offset(log(expected_2020)),
## family = "poisson", data = model_data_2020)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -78.13 -29.63 10.56 25.70 86.49
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.2336361 0.0075204 -31.07 <2e-16 ***
## prop_pop_3km 0.0048908 0.0001526 32.04 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 37822 on 23 degrees of freedom
## Residual deviance: 36775 on 22 degrees of freedom
## (3 observations deleted due to missingness)
## AIC: 37018
##
## Number of Fisher Scoring iterations: 5sjPlot::tab_model(model_3km_2020)| observed_2020 | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| (Intercept) | 0.79 | 0.78 – 0.80 | <0.001 |
| prop_pop_3km | 1.00 | 1.00 – 1.01 | <0.001 |
| Observations | 24 | ||
| R2 Nagelkerke | 1.000 | ||
Check how well the fitted values line up with the observations
The fitted values appear to line up particularly well with the observed data, suggesting that prop_pop_* i.e., proportion of catchment population living near water bodies can help us understand malaria risk in the catchment areas.
# Helper function to create scatter plots to see how well
# fitted values line up with observed malaria cases
plot.fitted.values <- function(fitted.values.df, model.df, title){
# Remove missing values from model data since
# model fitting deletes missing observations
model.df.complete <- model.df %>%
tidyr::drop_na() %>%
dplyr::rename_at(vars(starts_with("observed_")), ~ str_c("observed"))
# Plot fitted versus observed values
scatter.plot <- ggplot2::ggplot()+
ggplot2::geom_point(aes(fitted.values.df$fitted.values,
model.df.complete$observed))+
ggplot2::theme_classic()+
ggplot2::labs(y = "Observed values",
x = "Fitted values",
title = title)
return(scatter.plot)
}
# Invoking function
# 2017 -------------------------------------------------------------------------
fitted_1km_2017 <- plot.fitted.values(model_1km_2017, model_data_2017, "2017: 1km model")
fitted_2km_2017 <- plot.fitted.values(model_2km_2017, model_data_2017, "2017: 2km model")
fitted_3km_2017 <- plot.fitted.values(model_3km_2017, model_data_2017, "2017: 3km model")
# 2018 -------------------------------------------------------------------------
fitted_1km_2018 <- plot.fitted.values(model_1km_2018, model_data_2018, "2018: 1km model")
fitted_2km_2018 <- plot.fitted.values(model_2km_2018, model_data_2018, "2018: 2km model")
fitted_3km_2018 <- plot.fitted.values(model_3km_2018, model_data_2018, "2018: 3km model")
# 2019 -------------------------------------------------------------------------
fitted_1km_2019 <- plot.fitted.values(model_1km_2019, model_data_2019, "2019: 1km model")
fitted_2km_2019 <- plot.fitted.values(model_2km_2019, model_data_2019, "2019: 2km model")
fitted_3km_2019 <- plot.fitted.values(model_3km_2019, model_data_2020, "2019: 3km model")
# 2020 -------------------------------------------------------------------------
fitted_1km_2020 <- plot.fitted.values(model_1km_2020, model_data_2020, "2020: 1km model")
fitted_2km_2020 <- plot.fitted.values(model_2km_2020, model_data_2020, "2020: 2km model")
fitted_3km_2020 <- plot.fitted.values(model_3km_2020, model_data_2020, "2020: 3km model")
# Layout scatter plots ---------------------------------------------------------
cowplot::plot_grid(fitted_1km_2017, fitted_1km_2018, fitted_2km_2018, fitted_3km_2018,
fitted_1km_2019, fitted_2km_2019, fitted_1km_2020, ncol = 3, nrow = 4)
Fig 10. How well the percentage of catchment population living around water bodies explain observed malaria incidence
Test for residual spatial autocorrelation
Here we’ll use adjacency as criterion
# Find adjacent polygons,